package core;
import java.util.*;

/** 
 * @opt all
 * @hidden
 * @opt inferdeppackage
 * @navassoc - - - core.Component
 * @navassoc - - - core.Specie
 */
class UMLOptions {}

/**
 * A chemical system is a collection of equilibria along with the
 * parameters such as mass balance constraints and temperature that is
 * to be solved.
 *
 * @author Charlie Moore
 * @
 */
public class ChemicalSystem extends Observable {
    
    // Name of this particular System
    private String name;

    // A list of all of the equilibria included in this system
    private Vector<Equilibrium> equilibria = new Vector<Equilibrium>();

    // A list of all of the species included in this system
    private Vector<Specie> species = new Vector<Specie>();

    // A list of all of the components included in this system
    private Vector<Component> components = new Vector<Component>();

    // The temperature of the system in celsius
    private double temp;

    // The ionic strength of the system
    private double ionicStrength;

    // The number of iterations to allow the solve algorithm to converge
    // before giving up
    private int maxIterations;

    // If the error in the solve argorithm is less than this number, it
    // converged.
    private double errorTolerance;

    
    
    //=================================
    // Getters and Setters
    //---------------------------------
    
    public String getName() { return name; }
    public void   setName(String newName) { 
        name = newName;
        setChanged();
        notifyObservers();
    }

    public double getTemp() { return temp; }
    public void   setTemp(double newTemp) { 
        temp = newTemp;
        setChanged();
        notifyObservers();
    }

    public double  getIonicStrength() { return ionicStrength; }
    public void setIonicStrength(double newStrength) { 
        ionicStrength = newStrength;
        setChanged();
        notifyObservers();
    }

    public int  getMaxIterations() { return maxIterations; }
    public void setMaxIterations(int newMax) { 
        maxIterations = newMax;
        setChanged();
        notifyObservers();
    }

    public double getTolerance() { return errorTolerance; }
    public void setTolerance(double newTolerance) { 
        errorTolerance = newTolerance;
        setChanged();
        notifyObservers();
    }

    public Vector<Equilibrium> getEquilibria() { return equilibria; }
    public void addEquilibrium(Equilibrium eq) { 
        equilibria.add(eq);
        //Add the components from the equilibrium to the components list
        for(Component c : eq.getComponents())
        {
            addComponent(c);
        }
        for(Specie s : eq.getSpecies())
        {
            addSpecie(s);
        }
        setChanged();
        notifyObservers();
    }

    public Vector<Component> getComponents() { return components; }
    public void addComponent(Component comp) { 
        boolean found = false;
        for(Component temp : getComponents())
        {
            if(temp.getName().equalsIgnoreCase(comp.getName()))
            {
                found = true;
                break;
            }
        }
        if(!found)
        {
            components.add(comp);
            setChanged();
            notifyObservers();
        }
    }

    public Vector<Specie> getSpecies() { return species; }
    public void addSpecie(Specie spec) { 
        boolean found = false;
        for(Specie temp : getSpecies())
        {
            if(temp.getName().equalsIgnoreCase(spec.getName()))
            {
                found = true;
                break;
            }
        }
        if(!found)
        {
            species.add(spec);
            setChanged();
            notifyObservers();
        }
    }

    public void clearEquilibria() { 
        equilibria.clear();
        setChanged();
        notifyObservers();
    }

    public void clearSpecies() { 
        species.clear();
        setChanged();
        notifyObservers();
    }
    
    public void clearComponents() { 
        components.clear();
        setChanged();
        notifyObservers();
    }
    
    //=================================
    // Real methods
    //---------------------------------

    /** 
     * Constructor for a Chemical System.  The default temperature is 25 degrees
     * celsius.
     * 
     * @param systemName what this system will be called
     */
    public ChemicalSystem(String systemName) {
        setName(systemName);
        setTemp(25.0);
    }
    
    /** 
     * Constructor for a Chemical System.
     * 
     * @param systemName what this system will be called
     * @param systemTemp the temperature in celsius to use when
     * solving the system
     */    
    public ChemicalSystem(String systemName, double systemTemp) {
        setName(systemName);
        setTemp(systemTemp);
    }

    
    /** 
     * Solves an equilibrium problem within given mass balance
     * constraints for a given set of components.  The species
     * concentrations are calculated from free component activities
     * and log(K) values, then the jacobian matrix Jac is calculated
     * from these concentrations and the stoichiometric coefficients
     * in A.  New estimates of component activities G are then
     * calculated by solving the system Jac * G = TE where TE is the
     * total mass balance vector.  The process repeats until either a)
     * the mass balance errors are small or b) the maximum number of
     * iterations is exceeded.
     *
     * See also the MINEQL users guide by Westall and Morel for further 
     * illumination.
     * 
     * @param totalKnown the number of components for which the total
     *                   dissolved amount is known
     * @param numSpecies the number of species
     * @return boolean true if the system solved correctly, false otherwise
     */
    public boolean solve(int totalKnown, int numSpecies) {
        // Keeps track of the magnitude of summed contributions to mole balance
        double absSum = 0;
        double[][] jac = new double[totalKnown][totalKnown];
        double[] freeX = new double[totalKnown];
        double tempSum;
        int i = 0;
        int com = 0, spe = 0;

        while (i < this.maxIterations) { 
            // Calculate log concetrations of each species from the
            // sum of the log(K) of formation + the sum of log(C) of
            // each component * the stoichiometric coefficient
        }
        return true;
    }
    
    
    
    /**
     * Solves the simultaneous linear equations A * X = B, where A is the
     * coefficient matrix, B is the constant vector and X is the solution 
     * vector the matrix A is n x n, and the vectors are n units each.
     * This routine is from Pizer and Wallace, 1983, Chapter 6, and uses
     * Gaussian substitution with scaling and pivoting to reduce error
     *
     * @param A   the coefficient matrix
     * @param B   the constant vector
     * @param n   the length of the units
     * @return X   the solution vector
     * @throws Exception throws a generic exception if the matrix arithmetic 
     *                   has a problem.
     */
    private double[] matrixSolve(double[][] A, double[] X, double[] B, int n) throws Exception {
        double[] scale = new double[n];
        // Permutation Index
        int[] p = new int[n]; 
        // Largest scaled pivot
        double max = 0;
        // Row containing the max
        int maxRow = 0;
        // Magnitude of pivots
        double scaleValue = 0; 
        // Row Multiplier
        double aMult = 0;
        // P Holder
        int temp = 0;
        // Counters
        int i,j,k;
        
        
        // Initialize the permutation index and scaling vector
        for( i = 1; i <= n; i++ ) {
            p[i] = i;
            scale[i] = 0;
            
            for ( j = 1; j <= n; j++ ) {
                scale[i] = scale[i] + Math.abs(A[i][j]);
            }
        }
        
        // Triangularize with pivoting
        for( i = 1; i < n; i++ ) {
            max = 0;
            
            // Find the pivot
            for( j = 1; j <=n; j++ ) {
                scaleValue = Math.abs(A[p[j]][i] / scale[p[j]]);
                
                // Update max
                if( scaleValue > max) {
                    max = scaleValue;
                    maxRow = j;
                }
            }
            
            // Swap the rows
            temp = p[maxRow];
            p[maxRow] = p[i];
            p[i] = temp;
            
            // Eliminate subdiagonal elements on the ith column
            for( j = i+1; j <= n; j++) {
                if( A[p[i]][i] == 0) {
                    //TODO: Figure out the correct exception to throw
                    throw new Exception("MatrixSolve failed eliminating subdiagonal elements");
                }
                
                aMult = A[p[j]][i] / A[p[i]][i];
                for( k = i+1; k <=n; k++ ) {
                    A[p[j]][k] = A[p[j]][k] - aMult * A[p[i]][k];
                }
                
                B[p[j]] = B[p[j]] - aMult * B[p[i]];
            }
        }
        
        // Back Substitution
        for( i = n; i >= 1; i-- ) {
            X[i] = B[p[i]];
            
            for( j = i+1; j <=n; j++ ) {
                X[i] = X[i] - A[p[i]][j] * X[j];
            }
            
            if( A[p[i]][i] == 0 ) {
                throw new Exception("Matrixsolve failed in back substitution");
            }
            
            X[i] = X[i] / A[p[i]][i];
        }
        
        return X;
    }
    
    
    public void save_definition() {}

    public void load_definition() {}

}
